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1  INTRODUCTION 

This  is  the  final  report  on  the  research  project  conducted  at  the  CPAT  for  the  US  Airforce,  on 
Hall  thrusters,  called  stationary  plasma  thrusters  (SPTs)  in  the  following.  The  aim  of  the 
project  is  to  compare,  evaluate,  and  possibly  improve  the  following  two  SPT  models: 

1)  The  two-dimensional  hybrid  model  developed  at  the  MIT  between  1994  and  1998,  by 
Mike  Fife  and  Manuel  Martinez-Sanchez,  and  documented  in  [Fif98]. 

2)  The  two-dimensional  hybrid  model  developed  at  the  CPAT  from  2000  to  present,  by  the 
authors  of  this  report,  and  documented  in  [HagOl]  and  [Hag02]. 

Both  models  attempt  to  provide  a  complete  simulation  of  the  temporal  and  spatial  behavior  of 
the  SPT  discharge,  combining  a  particle  description  of  neutrals  and  ions  with  a  fluid 
description  of  electrons,  where  the  electric  field  is  obtained  from  assuming  quasi-neutrality. 

In  the  first  and  main  part  of  this  report  we  break  the  models  down  into  components,  which  we 
describe,  compare,  and  discuss  in  detail.  We  consider  the  following  topics: 

1)  Geometry  and  grid . 2 

2)  Magnetic  field  and  stream  function . 3 

3)  Neutral  atoms . 3 

4)  Sampling  of  neutral  atoms . 3 

5)  Accommodation  to  the  wall  temperature . 4 

6)  Ions . 5 

7)  Basic  electron  equations . 5 

8)  Anode  region . 6 

9)  Downstream  region . 7 

10)  Cross  field  electron  mobility . 7 

11)  Rates  of  ionization  and  collisional  energy-loss . 8 

12)  Wall-collision  theory  of  the  MIT  model . 9 

13)  Some  remarks  on  the  wall  theory  of  the  MIT  model .  11 

14)  Wall  effects  in  the  CPAT  model . 11 

15)  Wall  effects:  local  approach  vs  volumetric  approach .  12 

16)  Wall  effects:  MIT  model  vs  CPAT  model . 12 

17)  Numerical  solution  of  the  electron  equations . 13 

In  the  second  part  of  the  report  we  check  for  gross  numerical  errors  in  the  models,  by  slightly 
modifing  them  so  as  to  represent  similar  physics,  and  then  comparing  their  results. 
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2  DESCRIPTION  AND  COMPARISON  OF  MODEL  COMPONENTS 


2.1  Geometry  and  grid 

The  MIT  and  CPAT  models  represent  only  the  axial  and  radial  dimensions  of  the  SPT 
geometry  and  discharge;  azimuthal  symmetry  is  assumed.  The  computational  domain,  which 
comprises  both  the  discharge  channel  and  the  near  exterior  of  the  thruster,  is  schematically 
shown  in  figure  1 ;  the  full  discharge  simulation  is  canned  out  only  within  a  certain  region  that 
is  confined  by  physical  walls  and  magnetic  field  lines. 

The  two  models  use  different  computational  grids;  see  figure  2.  The  MIT  model  has  a  non- 
uniform  conformal  grid  that  may  be  adapted  to  arbitrarily  shaped  channel  walls.  The 
rectangular  grid  of  the  CPAT  model  does  not  offer  this  possibility,  but  on  the  other  hand  leads 
to  much  simpler,  and  thus  faster,  computation.  Inside  the  channel  both  grids  are  about  equally 
refined;  typically  a  grid  cell  measures  1  mm  x  1  mm. 


channel 


exterior 


Figure  1.  Schematic  picture  of  the  computational  domain. 


Figure  2.  Comparison  between  the  computational  grids  of  MIT  model  on  the  left-hand  side  and  the 
CPAT  model  on  the  right. 
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2.2  Magnetic  field  and  stream  function 

The  magnetic  field  is  assumed  to  be  entirely  determined  by  the  electromagnets  and  to  be 
unaffected  by  the  discharge,  so  that  it  can  be  treated  as  input  data.  Both  models  calculate,  to 
facilitate  the  solution  of  the  electron  transport  equations,  a  magnetic  stream  function  A  from 

M. =rB  and  —  =  -rB  til 

8x  '  dr  x’  1  1 

where  x  and  r  are  the  axial  and  radial  position  coordinates  and  Bx  and  B,  are  the  axial  and 
radial  components  of  the  magnetic  field.  The  A  defined  by  these  equations  is  constant  along 
magnetic  field  lines  (B  VT  =  0)  and  usually  increases  monotonically  from  anode  to  cathode. 

The  cross  field  gradient  of  any  quantity  Q  can  be  expressed  in  terms  of  A  as 

<2> 

An  additional  computational  grid  is  attached  to  the  /-coordinate  with  intervals  that  have 
approximately  the  same  size  as  the  cells  of  the  x-r- grid. 


2.3  Neutral  atoms 

In  both  models  the  density  of  neutral  xenon  atoms  in  the  thruster,  essential  to  find  the 
ionization  rate  and  the  electrical  conductivity  of  the  plasma,  is  obtained  from  a  Monte  Carlo 
simulation.  That  is,  the  individual  paths  of  a  large  number  of  neutrals  are  calculated,  where 
collisions  are  treated  with  random  numbers.  This  approach  is  realistic  but  takes  much 
computation  time  and  introduces  statistical  errors.  The  neutrals  are  introduced  in  the 
simulation  at  a  certain  injection  region  at  the  anode  and  are  followed  until  they  reach  the  right 
boundary  of  the  geometry.  Additional  neutrals  are  introduced  at  the  channel  walls  to  account 
for  wall-recombination  of  ions.  The  initial  neutral  velocity  distribution  is  randomly  sampled 
from  a  half-Maxwellian  distribution;  see  section  2.4.  Only  collisions  with  walls  are 
considered,  in  which  the  neutrals  are  isotropically  scattered. 

In  both  models  neutral  loss  by  ionization  is  implemented  by  the  following  procedure:  To  each 
simulated  neutral  a  certain  weight  w  is  attributed,  which  gradually  decreases  in  time  as 
w=w0exp(-nk,t),  (3) 

where  w0  is  the  initial  weight,  n  is  the  local  plasma  density,  k,  is  the  ionization  rate  coefficient 
dependent  on  the  local  electron  mean  energy,  and  t  is  the  time.  This  technique  leads  to  better 
statistics  than  simply  eliminating  neutrals  from  the  simulation  according  to  their  ionization 
probability,  especially  beyond  the  exhaust  where  the  neutral  density  may  be  quite  low. 


2.4  Sampling  of  neutral  atoms 


Confusion  may  easily  come  about  on  how  to  properly  determine  random  initial  velocities  for 
the  neutral  atoms  in  the  Monte-Carlo  model.  The  MIT  and  CPAT  models  use  different  but  in 
fact  equivalent  methods.  Here  we  document  these  methods. 

The  neutrals  are  introduced  at  the  anode  per  unit  time  and  per  unit  surface;  they  are  thus 
sampled  from  the  neutral  flux  crossing  the  anode  plane,  i.e.  from  the  neutral  flux  in  the  axial 
direction.  For  a  Maxwellian  velocity  distribution  the  axial  flux  in  the  Cartesian  velocity 
interval  dvxdv,,dv,  around  (vx,  vy,  vz)  is 

f 


&(c,vy,vj  dvvdvvdv,  cc  vxexp 


Mv 


2kT 


exp 


My, 


2kT 


exp 


Mv. 


2kT 


dvxdvydvz 


(4) 


where  k  is  the  Boltzmann  constant,  T  is  the  atom  temperature,  and  M  the  atom  mass. 
The  MIT  models  samples  directly  from  this  distribution  as 
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v-  = 

(5) 

Vy  =  -J|jln(^2)  COs(2/t/C) 

(6) 

A  =  ^HR2)  sin(2 /ri?3) 

(7) 

where  R \ _3  are  random  numbers  on  [0,1];  see  also  the  book  of  Bird  [Bir94].  On  average  vx  is 
larger  than  vv  and  vz  because  we  sample  from  an  axial  flux,  not  from  a  density. 

In  an  alternative  approach  the  CPAT  model  uses  the  distribution  in  spherical  coordinates 

mv~ 

2kT 


where  v  is  the  speed,  0  the  angle  with  the  axial  direction,  and  <j)  an  angle  in  the  plane 
perpendicular  to  the  axis.  Random  values  for  these  coordinates  now  follow  from 


hTrM-XrJ  -  * 

(9) 

COS0  =  yfll 7 

(10) 

<j)  -  2  xR3 

(ID 

where  i?i_3  are  once  again  three  random  numbers.  The  problem  here  is  that  the  equation  for  the 
speed  v  has  no  analytical  solution  and  requires  a  few  Newton  iterations.  On  the  other  hand, 
the  equations  for  6  and  (f)  are  extremely  useful  since  they  apply  not  only  to  a  Maxwellian 
distribution,  but  to  any  isotropic  distribution;  for  instance,  they  can  be  used  together  with  a 
fixed  value  for  the  speed  to  simulate  mono-energetic  particles. 

Analogous  formulas  are  being  used  for  neutrals  scattered  isotropically  at  the  walls  or  coming 
off  the  walls  from  recombination.  Note  that  in  this  case  the  radial  direction  takes  the  role 
played  in  the  above  formulas  by  the  axial  direction. 

Typically  the  MIT  model  uses  a  higher  initial  neutral  temperature  than  the  CPAT  model:  1000 
K,  against  500  K  for  the  CPAT  model.  The  MIT  value  is  based  on  measurements  of  the  anode 
temperature. 

2.5  Accommodation  to  the  wall  temperature 

The  MIT  model  assumes  that  the  neutrals  accommodate  to  the  wall  temperature  when 
colliding  with  the  channel  walls.  The  neutral  temperature  after  a  wall  collision  is  calculated  as 

^kTa=a^kTw+(l-a)^Mv2  (12) 

where  a  is  the  accommodation  coefficient,  Tw  is  the  wall  temperature,  and  Mv2/ 2  is  the  energy 
of  the  incident  neutral;  typically  a  =  0.8  and  Tw  =  900  K  constant  along  the  channel.  The 
neutral  velocity  after  the  collision  is  draw  randomly  from  a  half-Maxwellian  distribution  at 
temperature  Ta. 

The  CPAT  model  does  not  include  wall  accommodation.  Wall  collisions  only  change  the 
direction  of  the  neutral  velocity,  not  its  magnitude.  A  random  direction  is  found  from 
equations  (10)  and  (11). 

Note  that  even  for  a  =  0  the  MIT  model  differs  from  the  CPAT  model.  We  remark  that  the 
wall-accommodation  procedure  of  the  MIT  model  is  odd:  it  involves  treating  a  single 
particle’s  energy  as  if  it  were  a  mean  energy,  which  might  lead  to  undesirable  artifacts.  It 
would  be  better  to  treat  the  accommodation  as  follows:  First  find  a  random  speed  vw  from  a 
half-Maxwellian  distribution  at  the  wall  temperature;  this  can  be  done  either  directly  from 
equation  (9)  or  indirectly  from  three  components  (5)-(7).  Subsequently  the  accommodated 
neutral  speed  after  the  wall  collision  can  be  taken  as 


g(v,cos0,</>)  v^sint?  dvd^d^  oc  v3exp^ 


sint?cost?  dvddd</) 
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with  v  the  incident  speed,  and  the  direction  can  be  found  from  equations  (10)  and  (11),  as  in 
the  CPAT  model. 


2.6  Ions 

Like  the  neutrals,  the  ions  are  in  both  models  described  by  a  Monte  Carlo  simulation.  They 
are  introduced  in  the  simulation  at  positions  that  are  randomly  chosen  according  to  the 
ionization  rate  profile.  The  initial  ion  velocity  distribution  is  isotropic  and  Maxwellian  at  the 
neutral  gas  temperature.  The  ions  are  assumed  to  be  accelerated  by  the  electric  field  only,  i.e. 
to  be  insensitive  to  the  magnetic  field.  They  are  followed  until  they  reach  any  of  the 
boundaries  of  the  simulation  domain;  ions  striking  the  walls  are  thus  assumed  to  recombine  at 
the  surface.  Besides  the  ion  density,  the  ion  Monte  Carlo  simulation  yields  the  ion  flux  and 
the  ion  energy  distribution. 

The  MIT  model  distinguishes  between  singly  and  doubly  charge  ions,  whereas  the  CPAT 
model  assumes  all  ions  to  be  singly  charged.  Although  the  doubly  charged  ions  in  the  MIT 
model  may  make  up  for  about  1 0%  of  the  total  ion  flux,  the  double  ionization  seems  to  have 
little  influence  on  the  discharge  in  the  channel. 


2.7  Basic  electron  equations 


The  models  describe  the  electrons  by  a  set  of  fluid  equations,  which  we  document  in  this 
section.  We  use  the  formalism  of  the  CPAT  model;  although  the  equations  may  appear 
different  from  those  in  Mike  Fife’s  these,  we  assure  that  they  are  identical.  Differences 
between  the  models  are  discussed  in  later  sections. 

In  the  fluid  approach  the  behavior  of  the  electron  density,  flux,  and  mean  energy  is  described 
by  the  first  few  moments  of  the  Boltzmann  equation  (transport  equations);  this  incorporates 
many  assumptions  and  is  not  entirely  realistic.  In  view  of  the  high  plasma  density  in  SPTs  it 
is  assumed  that  the  electron  density  is  everywhere  equal  to  the  ion  density.  With  this 
assumption  it  becomes  impossible  to  obtain  the  electric  field  from  Poisson’s  equation. 
Instead,  knowing  the  electron  density,  we  use  the  electron  transport  equations  to  calculate  the 
electric  field. 

The  electron  transport  equations  are:  the  continuity  equation 

VTe  =  Nnkj  -  ^ — VT , ,  (14) 


(15) 


(16) 


the  momentum  equation,  which  we  approximate  by  the  drift-diffusion  equation 
Te=-juEn-^juV(ne), 
and  the  energy  equation 

+  |V-(rt,£-) - ^V  (/jnsVs) = -eETe  -Nn/c-nW  . 

In  these  equations  n  is  the  plasma  density,  Te  the  electron  flux,  s  the  electron  mean  energy,  N 
the  gas  density,  T,  the  ion  flux,  E  the  electric  field,  //  the  electron  mobility,  and  e  the 
elementary  charge.  The  last  two  terms  in  the  energy  equation  represent  energy  loss  by 
collisions  with  gas  particles  and  with  the  walls,  respectively,  where  k  and  W  are  effective 
energy  loss  coefficients  dependent  on  s.  The  equations  (15)  and  (16)  assume  the  electron 
distribution  to  be  Maxwellian  and  predominantly  isotropic;  the  same  assumption  is  used  to 
obtain  the  collision  coefficients  k,  and  k  from  cross  section  data. 

Due  to  the  magnetic  field  the  mobility  //  is  not  a  simple  scalar:  its  value  is  much  larger  for 
electron  transport  along  magnetic  field  lines  than  for  transport  across  them.  From  current 
conservation  however  it  is  clear  that  the  electron  flux  cannot  be  much  larger  along  magnetic 
field  lines  than  across  them.  This  implies  that  along  the  field  lines  the  two  terms  of  equation 
(15)  should  virtually  cancel  each  other.  Taking  into  account  similar  considerations  for  the 
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electron  energy  flux,  one  can  derive  that  the  electron  mean  energy  must  be  constant  along 
magnetic  field  lines  and  that  the  electric  potential  V  behaves  as 


n(x,r ) 


V{x,r)=V\X)+-^4^) In-  no 


(17) 


where  V  is  a  function  and  «0  is  a  reference  density.  While  V  and  n  vary  all  over  space,  V  and 
s  depend  only  on  the  stream  function  A.  Note  that  by  using  this  equation  we  lose  to  possibility 
to  calculate  the  electron  flux  along  magnetic  field  lines  from  the  drift-diffusion  equation.  For 
the  cross  field  electron  flux  on  the  other  hand  we  find 


r  ~rBu  n^-——rBu  ^ n£  ^ 
\eA-rB^ndx  3g rB /U  ^ 


=  rB/i1n@j!fr  +  ^rBn±n I  ln-^ — 1  l-g^- 


dA  3e 


ds 


dA 


where  /j±  is  now  the  cross  field  mobility,  which  is  further  discussed  in  Section  lll.A. 
Let  us  now  define  the  following  (surface!)  integrals  along  magnetic  field  lines 
t  i  I ,  i  ds 

=||rS//1«dv 


c2 


c3: 


In— — 1 
«o 


(18) 

(19) 

(20) 

(21) 


and  (volume)  integrals  between  consecutive  field  lines 

c4=|||«dv  (22) 

c5=[f|/V«dv’  (23) 

c^JJ-^dv,  (24) 

where  d.v  and  dv  are  surface  and  volume  elements. 

Using  these  integrals,  the  continuity  and  momentum  equations  can  be  replaced  by  the 
following  one-dimensional  equation  for  current  conservation: 

(25) 


jjf,,d.s  =c2 


dV_  .2.  ds_r  _1  j 
QA  3  e  3  dA  1  e  ’ 


where  /  is  the  discharge  current.  It  is  assumed  that  no  current  escapes  to  the  walls. 
In  a  similar  way,  the  energy  equation  can  be  written  as 

5i  _\a  ^ 


^  4  A  k  ^  +  -|(ci,A-+l/2  “  ^A+l/2-f  pk+l/2 


dt 

10,,  „  ds 

■  9eC2MU2£k+y  2dA 
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,10  ,.  ds 

+  9eC2,k-V2£k- V2qa 


--c6,k  -c5'kKk-c*kWk 


(26) 


where  k+ 1/2  and  k- 1/2  refer  to  two  field  lines,  and  k  to  the  interval  between  them. 

From  the  equations  (25)  and  (26)  we  calculate  s  and  V  as  a  function  of  /..  Subsequently  the 
spatial  profile  of  the  electric  potential  is  found  from  equation  (17).  The  current  in  the 
equations  (25)  and  (26)  is  chosen  such,  that  a  specified  voltage  drop  results  between  anode 
and  cathode: 


K-V=- 


‘  4^-d A  +  ^-sa  In—  -  ^-sc  In— , 
i  oA  3e  »n  3e 


(27) 


n0  ie  n0 

where  the  labels  a  and  c  refer  to  anode  and  cathode,  and  the  naiC  are  spatially  averaged  over 
the  electrode  lines. 


2.8  Anode  region 

The  physics  of  the  region  in  front  of  the  anode  is  not  well  described  by  the  above  electron 
equations,  which  are  solved  only  right  of  anode  line,  i.e.  the  first  magnetic  field  line  from  the 
left  that  does  not  intercept  the  anode.  The  electron  mean  energy  and  electric  potential  on  the 
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anode  line  and  in  the  region  left  of  it,  where  the  magnetic  field  intercept  the  anode,  result 
from  boundary  conditions  and  assumptions. 

The  MIT  model  imposes  a  zero  energy  gradient  on  the  anode  line;  left  of  it  the  MIT  model 
assumes  the  electrons  to  be  fully  accommodated  to  the  anode  temperature  (1000  K).  This 
leads  to  an  unnatural  discontinuity  of  the  electron  energy  on  the  anode  line;  the  function  V  is 
ajusted  so  as  to  minimize  the  discontinuity  in  the  potential  profile. 

The  CP  AT  model  simply  fixes  the  electron  energy  to  a  small  value,  typically  on  the  order  of  1 
eV,  not  only  on  the  anode  line  but  also  in  the  region  left  of  it. 

We  remark  that  in  either  model  too  high  electron  mean  energy  at  the  anode  (>  7  eV)  leads  to 
excessive  bulk  ionization  there,  even  a  global  ionization  maximum,  which  seems  hardly 
realistic.  In  the  MIT  model  this  happens  if  the  magnetic  field  is  too  strong  at  the  anode;  in  the 
CP  AT  model  the  usual  boundary  condition  prevents  it  from  happening.  From  a  physical  point 
of  view,  one  might  imagine  that  collisions  with  anode  cool  the  electrons. 


In  contrast  to  the  CPAT  model,  the  MIT  model  accounts  for  the  voltage  drop  across  the  anode 
sheath  when  calculating  the  potential  distribution  from  the  electron  equations.  The  anode 
sheath  voltage  is  estimated  from  the  classical  expression 


where  kTe  =  (3 12)  e,  me  is  the  electron  mass,  m,  is  the  ion  mass,  A  is  the  surface  area  of  the 
anode,  vth  =  (8 kTJjm^'1,  and  e  =cxp(  1 );  this  is  only  valid  for  ion-attracting  sheaths,  i.e.  Vs 
positive. 


2.9  Downstream  region 

Both  the  MIT  and  CPAT  models  fix  the  electron  mean  energy  to  a  few  eVs  on  the  cathode 
line,  i.e.  the  magnetic  field  line  that  intercepts  the  anode.  For  the  region  downstream  of  the 
cathode  line  there  is  a  crucial  difference  between  the  models,  which  we  point  out  below. 

The  MIT  model  does  not  solve  the  electron  equations  in  the  downstream  region.  Instead,  it 
defines  an  effective  ground  point,  located  at  the  utter  right  egde  of  the  grid,  where  it  sets  s  and 

V  to  experimental  values  for  the  far-field;  between  the  cathode  and  the  effective  ground  e  and 

V  are  then  linearly  interpolated  on  the  /.-grid.  Obviously  this  technique  does  not  describe  the 
physics  of  the  downstream  region;  it  should  be  regarded  as  an  attempt  to  compress  the  entire 
downstream  region  into  the  computation  domain,  so  as  to  obtain  a  more  accurate  prediction 
for  the  eventual  ion  exhaust-velocity. 

The  CPAT  model  on  the  other  hand  solves  the  full  set  of  electron  equations  in  the 
downstream  region.  Given  that  for  each  ion  leaving  the  thruster  an  electron  must  leave  it  too, 
the  total  current  in  the  current-conservation  equation  (25)  is  set  to  zero  here.  On  the  right  edge 
of  the  calculation  domain  the  electron  mean  energy  is  fixed  to  a  small  value.  Typically  the 
CPAT  model  predicts  a  rise  in  the  electric  potential  on  the  order  of  1 0  V  beyond  the  cathode 
line;  further  downstream  the  potential  decreases  again;  this  finding  is  supported  by  the 
experimental  data  in  Fife’s  thesis  [Fif98],  Accurate  prediction  of  the  final  ion  exhaust- 
velocity  requires  a  large  computation  domain. 

2.10  Cross  field  electron  mobility 

Cross  field  electron  mobility  is  the  main  parameter  controlling  the  electric  potential 
distribution  in  the  SPT  and  has  therefore  an  enormous  influence  on  the  simulation  results. 
[HagOl]  Unfortunately  it  is  not  well  known;  the  MIT  and  CPAT  models  use  rather  different 
assumptions  for  it. 

The  classical  expression  for  /u±  is  given  by 
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(29) 


..  _  evjme  _rnevm 
v.'tW  eB2  ■ 

where  me  the  electron  mass  and  vm  the  momentum-transfer  frequency  of  electron-particle 
collisions.  The  MIT  model  assumes  a  constant  momentum-transfer  cross-section  of  3xl0‘19 
m2  whereas  the  CPAT  takes  the  momentum-transfer  frequency  to  be  constant  at  2.5xl0'13  m3 
s'1;  see  figure  3;  the  CPAT  model  is  closer  to  reality. 

It  is  a  known  fact  that  the  classical  mobility  is  too  small  to  be  realistic  for  the  electron 
transport  in  SPTs,  especially  near  and  beyond  the  exhaust  where  the  gas  density  is  very  low. 
Both  models  therefore  add  to  the  classical  mobility  an  anomalous  Bohm  mobility 

^=r h  <30) 

where  K  is  a  constant  fitting  parameter.  Anomalous  Bohm  mobility  has  been  measured  in 
various  magnetized  plasma’s  and  is  usually  physically  interpreted  as  the  result  of  plasma-field 
fluctuations. 

The  MIT  model  applies  Bohm  mobility  with  a  factor  K  =  0.15  in  the  entire  computation 
domain.  Doing  so,  the  model  predicts  the  acceleration  zone  to  be  entirely  located  outside  the 
channel;  see  CASE  I  in  Fife’s  thesis  [Fif98].  Since  this  is  clearly  not  realistic,  the  MIT  model 
moves  the  cathode  (now  called  ‘effective  cathode’)  very  close  to  the  channel  exhaust;  this 
pushes  the  acceleration  zone  in;  see  CASE  II  in  Fife’s  thesis.  Note  that  in  the  MIT  model  the 
electric  potential  profile  downstream  of  the  cathode  results  from  interpolation  instead  of  from 
electron  equations;  it  follows  that  the  exterior  is  not  really  modeled  at  all. 

The  CPAT  model  argues  that  to  have  the  acceleration  zone  inside  the  channel,  the  electron 
mobility  must  be  much  larger  outside  than  inside.  It  therefore  uses  the  Bohm  mobility  only 
outside,  where  K  =  0.2-1;  the  acceleration  zone  is  indeed  found  inside.  It  has  been  suggested 
in  the  litterature  that  the  plasma  flow  be  more  stable  inside  due  to  the  gradient  of  the  magnetic 
field  strength  [Zhu99],  which  supports  the  CPAT  treatment  of  Bohm  mobility.  The 
discontinuity  in  the  CPAT  mobility  at  the  channel  exhaust  seems  to  have  no  serious 
consequences  for  the  simulation  results. 

In  addition  to  classical  and  Bohm  mobility,  electron-wall  collisions  contribute  to  cross-field 
electron  transport.  This  is  described  in  later  sections.  The  various  mechanisms  of  electron 
transport  are  compared  in  figure  6. 


0  20  40  60  80 

electron  mean  energy  (eV) 

Figure  3.  Comparison  of  vJN,  where  N  is  the  neutral  density,  assumed  by  the  two  models. 


2.11  Rates  of  ionization  and  collisional  energy-loss 

The  electron  equations  involve  the  ionization  rate  coefficient  k,  and  the  collisional  energy-loss 
coefficent  k,  both  as  a  function  of  electron  mean  energy.  These  coefficients  result  from 
integrating  collision  cross-sections  over  a  Maxwellian  electron  energy  distribution  function. 
The  MIT  model  uses  theoretical  data  from  Drawin  and  Dugan  [Dug],  the  CPAT  model 
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measured  data  from  Puech  and  Mizzi  [Pue91].  The  resulting  coefficients  are  somewhat 
different;  figure  4  gives  a  comparison.  The  MIT  model  has  a  smaller  ionization  coefficient  but 
a  larger  energy-loss  coefficient.  We  remark  that  it  seems  odd  that  according  to  the  MIT  data, 
high  energy  electrons  undergo  an  effective  energy  loss  of  about  twice  the  ionization  threshold 
per  ionization  event. 


Figure  4.  Comparison  of  the  collision  coefficients  k,  and  /fused  by  both  models. 


2.12  Wall-collision  theory  of  the  MIT  model 


The  MIT  model  attempts  to  give  a  quantitative  description  of  the  effects  of  electron-wall 
collisions,  taking  into  account  electron-impact  secondary  emission.  Here  we  reconstruct  the 
MIT  wall-theory  from  the  information  in  Mike  Fife’s  thesis  [Fif98]. 

Experimental  data  on  secondary  emission  are  fitted  as 

sJ^-=AEb  (31) 

where  F(,,  and  Tep  are  secondary-electron  and  (mono-energetic)  primary- electron  fluxes,  and 
E  is  the  (mono-energetic)  primary-electron  impact  energy.  The  fitting  parameters  A  and  B  are 
typically  set  to  0.10-0.15  and  0.6,  respectively.  Note  that  here  the  secondary-electron  flux 
represents  all  electrons  coming  off  the  wall,  including  the  back-scattered  primaries.  By 
integration  over  a  Maxwellian  distribution  for  the  primary-electron  energy  E  the  expression  is 
found  for  the  effective  secondary-emission  coefficient 

S  =  AT{2+B){kTelef ,  (32) 

where  T  is  Euler’s  gamma-function  and  Te  is  the  temperature  of  the  primaries;  in  the 
following  it  is  assumed  that  the  primaries  at  the  wall  and  electrons  in  the  plasma  have  the 
same  temperature  Te;  this  is  valid  for  Boltzmann  equilibrium  in  the  sheath. 

Subsequently  the  wall  voltage  Vw  with  respect  to  the  plasma  is  derived  from  imposing  zero 
net  current  to  the  surface.  For  a  normal  ion- attracting  sheath  the  following  balance  equation  is 
considered  at  the  sheath  edge: 

,i .... 


kTp 


-nvR 


(33) 


/ 

\  1/2 


where  vth  =  (SkTJ  7me)lu  is  the  electron  thermal  velocity,  vB  =  ( kTJm ,)m  is  the  ion  Bohm 
velocity,  and  Vp=kTe/e  is  the  presheath  voltage  as  predicted  by  the  standard  Bohm  theory.  This 
yields 

kT 


K,  = 


-In 


M 


em, 
2  mn„ 


(34) 


where  e  =cxp(  1 ).  If  5  approaches  1  the  sheath  voltage  goes  positive  and  equation  (33)  is  no 
longer  valid.  Instead,  the  wall  voltage  for  an  ion-repelling  sheath  found  from  balancing  the 
primary-  and  secondary-electron  fluxes,  neglecting  the  ion  flux: 
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j-/7v„,exp| 


eVw 

kZ 


e,s  J 


that  is 

K  = 


kz 


ln<7. 


(35) 


(36) 


Here  Tes  is  the  temperature  the  secondary  electrons,  which  is  typically  assumed  to  be  on  the 
order  of  the  wall  temperature  (0. 1  eV).  Note  that  the  transition  between  the  two  regimes  is 
uniquely  characterized  by  the  break-point  temperature 

'he — 


T.h  = 


1 7imejemi 


(37) 


for  which  <7 is  close  to  1.  Typically  Teb  =15-30  eV. 

In  the  above  sheath  model  the  primary-electron  flux  to  the  wall  is 


r  =• 

e,p 


^nvlhex 


eV 

^  r  w 

kZ 


Z<Z 


e,b 


(38) 


1e^1e,b 

where  the  wall  voltage  for  Ze<Zeb  is  given  by  equation  (34). 

The  MIT  model  goes  on  arguing  that  the  secondary  emission  leads  to  a  net  cross-field 
electron  transport,  because  the  eventual  orbits  of  the  secondary-electrons  are,  on  average, 
shifted  towards  the  anode  by  one  Larmor  radius  with  respect  to  the  primary-electron  orbits. 
The  resulting  wall-current  is  shown  to  satisfy  the  equation 
2  7wmeE 


7>7I, 


I^eTepS- 


eBl 


where  r  is  the  radius  of  the  wall.  Substituting  some  of  the  above  equations  one  finds 

Ze<Zeb 


EM 


jkZ  5  2m-Eme 

\emi  1-5 

B2 

1  kZe  ,,  2mZme 

\27nne 

B2 

(39) 


(40) 


z>z, 


The  MIT  model  also  predicts  the  electron- energy  flux  to  the  wall.  For  an  ion-attracting  sheath 
the  primary  electrons  lost  at  the  wall  had,  before  they  enterred  the  sheath,  a  mean  energy 
equal  to  2 kZe-eVw  (Vw  negative!);  the  emitted  secondaries,  once  they  reach  the  plasma,  have 
been  accerelated  to  a  mean  energy  of  2kTes-e  Vw.  The  net  energy  flux  is  therefore 
qw={2kZ-eVw)re,P  ~  {2kTetS-eVw)ST^  Te<TeJb .  (41) 

For  an  ion-repelling  sheath  on  the  other  hand,  the  primary  electrons  are  ordinary  plasma- 
electrons  with  mean  energy  2kZe  and  the  secondary  electrons  that  make  it  to  the  plasma  arrive 
there  with  mean  energy  2 kZe  s  so  that  the  wall  voltage  does  not  interfere  in  the  expression 

qw=2kZtre,p-2kZe  Ze,P  Te>Ze^b . 

Using  the  above  equations  one  ends  up  with 


Mi 


em , 
2/Z777„ 


71  <71 


e,b 


J) 


(42) 


z>z 


e,b 


The  wall  current  (40)  and  the  energy  loss  flux  (42)  are  explicitly  taken  into  account  in  the 
current  conservation  and  energy  equations,  respectively.  For  instance,  the  current 
conservation  equation  of  the  MIT  model  is 

Jjr.xdy  -Mi  ,  (43) 

J  *  c  C 

where  Iw\  and  /„,2  are  the  wall  currents  on  the  inner  and  outer  walls. 
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2.13  Some  remarks  on  the  wall  theory  of  the  MIT  model 

1)  The  wall  theory  depends  strongly  on  the  data  for  the  secondary-emission  coefficient. 
What  data  are  being  used?  Is  the  fitting  formula  really  appropriate?  The  data  referred  to  in 
Mike  Fife’s  thesis,  viz.  measurements  published  by  Bugeat  et  al.  at  the  1995  IEPC 
[Bug95],  look  rather  limited.  In  the  few  other  data  that  we  found  ourselves  the  mean 
energy  of  the  secondary  electrons  (including  reflected  primaries)  is  a  significant  part  of 
that  of  the  primaries,  if  not  close  to  it,  i.e.  much  larger  than  Te  s  ~  0. 1  eV  as  assumed  in  the 
MIT  model. 

2)  The  wall  theory  relies  heavily  on  the  assumption  of  a  Maxwellian  electron-energy 
distribution.  This  is  difficult  to  justify;  the  wall-collisions  themselves  can  be  expected  to 
completely  alter  the  electron  energy  distribution,  in  particular  to  deplete  its  tail.  As  a 
consequence  the  theory  might  strongly  overestimate  the  effects  of  wall-collisions.  In  fact, 
deviations  from  Maxwellian  do  not  only  affect  the  wall  theory  but  the  entire  electron 
model. 

3)  The  wall  voltage  of  an  ion-attracting  sheath  includes  the  presheath  voltage,  which  implies 
that  the  plasma  density  used  to  calculate  the  primary-electron  flux  to  the  wall  (38)  is  to  be 
taken  outside  the  presheath.  But  where  is  that?  Seeing  that  the  ions  are  virtually 
collisionless,  the  presheath  will  stretch  out  to  the  center  of  the  plasma;  this  wall  theory  is 
therefore  not  local  but  involves  the  entire  plasma  volume. 

4)  The  ion  model  predicts  the  ion  flux  to  the  wall,  albeit  with  strong  statistical  noise.  This 
information  could  be  used  in  the  wall  model  to  circumvent  the  Bohm  theory  and  obtain 
more  consistency.  Equation  (33)  could  then  be  replaced  by 

(l- 8 )^nvth  ex^j^j=T, ,  (44) 

and  so  on.  The  primary-electron  flux  for  an  ion-attracting  sheath  would  simply  be 
r^=T=jr'-  (45) 


2.14  Wall  effects  in  the  CPAT  model 


The  CPAT  model  does  not  venture  into  a  quantitative  description  of  electron-wall  interactions 
and  limits  itself  to  being  qualitative. 

Cross-field  electron  transport  due  to  wall  collisions  is  taken  into  account  by  an  additional 
contribution  to  the  cross-field  electron  mobility.  The  physical  interpretation  of  this  approach 
is  that,  except  for  the  position  were  they  occur,  there  is  no  essential  difference  between 
electron-wall  collisions  and  electron-particle  collisions.  The  contribution  to  the  mobility  from 
the  wall  collisions  is 


Mw=~ 


eB2 


(46) 


where  vw  is  the  momentum  transfer  frequency  of  the  wall  collisions,  which,  for  the  sake  of 
simplicity,  is  kept  constant  along  the  channel: 

vm  =al07  (s'1).  (47) 


The  fitting  parameter  a  is  typically  chosen  in  the  range  0.1 -0.5. 

Energy  loss  due  to  electron-wall  collisions  is  accounted  for  via  an  empirical  volumetric 
energy  loss  coefficient  (energy  loss  per  second  per  electron) 

W  =  a£  107  f  exp(— (s'1),  (48) 

where  as  and  U  are  once  again  constant  fitting  parameters,  typically  0. 1-0.5  and  20  eV, 
respectively.  This  does  not  result  from  physical  derivation,  but  is  simply  a  convenient 
expression  stating  that  the  energy  wall-loss  increases  monotonically,  becoming  important 
beyond  the  energy  U. 
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2.15  Wall  effects:  local  approach  vs  volumetric  approach 


The  MIT  model  uses  an  explicit  equation  for  the  wall  current  whereas  the  CPAT  includes  it 
volumetrically  via  the  electron  mobility.  Although  at  first  sight  these  approaches  may  look 
very  different,  they  are  not. 

For  demonstration  we  translate  the  MIT  expression  for  the  wall  current  into  an  equivalent 
mobilty  contribution.  Consider  a  cross  section  through  the  channel,  perpendicular  to  the  axis. 
Assume  that  the  magnetic  and  electric  fields,  the  plasma  density,  and  the  electron  temperature 
are  all  constant  over  this  cross  section.  The  wall  current  through  the  cross  section  in  terms  of  a 
wall  mobility  is  then 


e7i{r2  -i\)/uwEn=en{r2  ,  (49) 

where  r\.2  are  the  radii  of  the  walls.  Equating  this  to  the  sum  of  the  wall  currents  (40)  on  the 
inner  and  outer  walls  gives 


2  Is  5 
r2-/|  \3entj  1  -8 

2 
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z<\kTe,b 

s>^kTeb 


(50) 


where  we  substituted  kTe  =  (3/2)  £ 

Similarly,  the  MIT  energy  flux  to  the  wall  can  be  written  in  the  form  of  the  volumetric  energy 
loss  coefficient  W.  Integrating  the  energy  loss  over  a  slice  perpendicular  to  the  axis  and  with  a 
thickness  d  one  obtains  the  equation 

n(r22 -r^d  nW  -  2n{}\+r2)d  qw,  (51) 

which  yields 
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(52) 


Note  however  that  there  are  some  minor  differences  between  the  explicit  approach  of  the  MIT 
model  and  the  equivalent  volumetric  approach: 


1)  In  the  MIT  model  the  wall  effects  are  calculated  from  the  local  fields  and  plasma  density 
at  the  walls,  whereas  the  volumetric  treatment  averages  over  the  volume.  It  is  not  clear 
that  locality  is  an  advantage  here,  given  the  fact  that  the  presheath  extends  far  into  the 
plasma;  see  also  remark  3)  of  section  2.13. 

2)  The  MIT  wall-current  (40)  is  directly  proportional  to  the  electric  field,  whereas  with  the 
volumetric  treatment  wall  collisions  also  lead  to  a  diffusion  flux.  The  latter  does  not  seem 
less  realistic  than  the  first. 


2.16  Wall  effects:  MIT  model  vs  CPAT  model 

Using  the  volumetric  equivalents  presented  in  the  previous  section,  one  can  directly  compare 
the  MIT  wall-theory  with  the  empirical  CPAT  formulas.  This  is  done  in  figure  5  for  typical 
values  of  the  various  fitting  parameters. 

In  order  to  give  an  idea  of  how  the  wall-effect  relate  to  the  rest  of  the  model,  figure  6 
compares  them  with  competing  volume  effects,  discussed  in  previous  sections;  each  of  the 
plots  in  this  figure  is  based  on  the  same  data  for  the  magnetic  field,  neutral  density,  and 
electron  mean  energy,  obtained  with  a  standard  simulation  that  is  further  described  in  the 
chapter  on  code-to-code  validation. 
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Figure  5.  Comparison  between  the  wall-collision  frequencies  and  energy-loss  coefficients  of  the  MIT 
and  CP  AT  models.  On  the  left  the  functions  (47)  and  (50)  are  plotted,  on  the  right  the  functions  (48) 
and  (52).  The  fitting  parameters  are  set  to  standard  values:  A  =  0.108,  B  =  0.576,  es  =  0.1  eV  for  the 
MIT  model  and  a=  as=  0.2,  U  =  20  eV  for  the  CPAT  model. 


axial  position  (m) 


Figure  6.  Comparison  of  the  wall  effects  in  the  MIT  and  CPAT  models  with  competing  volume  effects. 
The  plots  show  axial  profiles  in  the  center  of  the  channel  for  a  typical  simulation,  which  is  described  in 
more  detail  in  chapter  3.  The  anode  cathode  are  located  on  the  left  and  on  the  right,  respectively.  The 
Bohm  parameter  K  =  0.15,  the  wall  fitting  parameters  are  set  to  the  same  standard  values  as  with  the 
previous  figure. 


2.17  Numerical  solution  of  the  electron  equations 

The  solution  of  the  electron  equations  is  tough  numerical  task  that  is  handled  very  differently 
by  the  two  models. 

The  MIT  model  solves  the  electron  energy  equation  by  a  modified  Forward  Time  Centered 
Space  method:  the  energy  is  updated  fully  explicitly  and  all  differentials  are  approximated  by 
central  differences  on  the  /.-grid.  This  technique  necessitates  the  use  of  a  very  small  time  step, 
typically  5x1 0'11  s,  so  that  many  iterations  must  be  performed  to  simulate  a  significant  period 
of  SPT  operation.  We  remark  that  central  differences  are  not  appropriate  for  the  first-order 
terms  in  space  if  they  dominate  the  second-order  terms;  in  practice  however  this  must  not  be 
the  case  since  the  MIT  method  turns  out  to  work  reasonably  well. 

The  CPAT  model  evaluates  the  spatial  differences  fully  implicitly,  using  a  very  robust 
exponential  discretization  scheme  [Sch69].  The  energy  loss  terms  are  evaluated  implicitly  by 
a  Newton  iteration,  i.e.  the  terms  are  linearized  in  energy  around  the  value  of  the  previous 
iteration.  [HagOO]  The  implicit  technique  avoids  strong  time  step  restictions;  time  steps  of  10"8 
s  or  more  can  be  used  without  problem.  As  a  result,  the  computation  time  spent  on  solving  the 
electron  equation  is  negligible  in  the  CPAT  model. 
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3  CODE-TO-CODE  VALIDATION 


In  order  to  check  for  gross  numerical  errors  in  the  models  and  thus  give  more  confidence  in 
their  results,  we  have  directly  compared  them  on  the  simulation  of  the  SPT-70.  For  this 
purpose  the  models  were  slightly  modified  so  as  to  represent  as  much  as  possible  the  same 
physical  equations.  The  modifications  included  the  following: 

1)  Double  ionization  was  turned  off  in  the  MIT  model. 

2)  The  MIT  wall-theory  was  included  in  the  CPAT  model  via  the  equivalent  wall-mobility 
and  energy  wall-loss  coefficient. 

3)  The  MIT  data  for  the  ionization  coefficient,  the  energy-loss  coefficient,  and  classical 
electron  mobility,  were  implented  into  the  CPAT  model. 

4)  Bohm  mobility  inside  the  channel  was  added  to  the  CPAT  model. 

5)  The  cathode  in  the  CPAT  model  was  placed  close  to  the  exhaust,  as  in  the  MIT  model. 

6)  The  potential  drop  across  the  anode  sheath  as  calculated  by  the  MIT  model  was  explicitly 
added  to  the  applied  voltage  in  the  CPAT  model. 

In  spite  of  these  changes  some  differences  were  still  present: 

1)  Since  the  CPAT  treatment  of  the  region  downstream  of  the  cathode  is  on  no  account 
compatible  with  the  MIT  assumptions  on  the  electron  mobility,  it  was  decided  to 
completely  turn  off  the  CPAT  simulation  of  this  region;  the  MIT  model  still  used  its  usual 
interpolation  technique.  Note  that  the  downstream  region  is  not  expected  to  have  much 
influence  on  the  discharge  inside  the  channel. 

2)  The  boundary  conditions  on  the  anode  were  left  as  is:  In  the  region  left  of  the  anode  line 
the  MIT  model  assumes  the  electron  mean  energy  to  be  0.1  eV,  whereas  the  CPAT  model 
uses  the  same  value  as  on  the  anode  line,  about  6  eV  in  this  case. 

3)  The  volumetric  CPAT  treatment  of  wall-effects  is  slightly  different  from  the  local  MIT 
treatment;  see  section  2.15. 

Standard  input  data  and  conditions  were  used:  magnetic  field  as  in  Fife’s  thesis,  applied 
voltage  300  V,  debit  2.34  mg/s,  Bohm  parameter  K  =  0.15,  wall  parameters  A  =  0.108,  B  = 
0.576,  £s  =  0. 1  eV. 

The  two  codes  turned  out  to  yield  very  similar  simulation  results  indeed,  in  spite  of  their  very 
different  numerical  methods;  see  figures  7  and  8.  Not  only  the  time-averaged  spatial  results 
are  close  together,  but  also  the  temperoral  behavior  and  calculated  oscillations.  The  minor 
differences  in  the  results  can  easily  be  attributed  to  the  above-mentionned  differences  that 
remained  in  the  physical  equations,  and  do  not  give  rise  to  doubts  about  the  numerical 
implementation  of  the  models. 


Figure  7.  Comparison  of  the  discharge  currents  calculated  by  the  MIT  and  CPAT  models,  modified  to 
represent  similar  physics,  for  a  standard  simulation  of  the  SPT-70. 
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Figure  8.  Direct  comparison  of  the  MIT  and  CPAT  models,  modified  to  represent  similar  physics,  on  a 
standard  simulation  of  the  SPT-70.  The  plots  show  time-averaged  two-dimensional  profiles  of  the 
neutral  density,  plasma  density,  ionization  rate,  electron  mean  energy,  and  electric  potential, 
respectively;  on  the  left-hand  side  the  MIT  model,  right  the  CPAT  model. 
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4  CONCLUSIONS 


There  are  two  major  differences  in  the  physics  of  the  MIT  and  CPAT  models: 

1)  The  MIT  model  estimates  the  effects  of  electron-wall  collisions  from  a  sheath  theory, 
taking  into  account  electron-impact  secondary-electron  emission,  whereas  the  CPAT 
model  uses  simple  empirical  formulas.  Although  the  MIT  model  approach  seems  more 
advanced,  it  is  based  on  many  questionnable  assumptions. 

2)  The  MIT  model  uses  the  same  Bohm  conductivity  everywhere,  while  CPAT  model  uses 
Bohm  conductivity  only  outside.  With  the  MIT  approach  the  cathode  must  be  placed  near 
the  exhaust  to  obtain  an  acceleration  zone  inside  the  channel;  the  region  beyond  the 
cathode  cannot  be  simulated  and  is  described  by  interpolation  to  a  ground  point.  The 
CPAT  approach  automatically  locates  the  acceleration  zone  inside  the  channel  and  allows 
the  full  simulation  of  the  downstream  region. 

These  differences  mark  the  points  were  the  physics  is  the  least  understood  and  the 
assumptions  are  the  most  doubtful.  Without  new  experimental  evidence,  significant 
improvement  of  the  models  seems  unlikely. 

When  modified  to  represent  similar  physics,  the  models  yield  very  similar  results.  Not  only 
the  time-averaged  spatial  results  are  close  together,  but  also  the  temperoral  behavior  and 
calculated  oscillations.  This  is  strong  evidence  that  there  are  no  gross  numerical  errors  in 
either  of  the  codes. 
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